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Abstract: Long-term trends in photosynthetic capacity measured with the satellite-derived 
Normalized Difference Vegetation Index (NDVI) are usually associated with climate 
change. Human impacts on the global land surface are typically not accounted for. Here, 
we provide the first global analysis quantifying the effect of the earth’s human footprint on 
NDVI trends. Globally, more than 20% of the variability in NDVI trends was explained by 
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anthropogenic factors such as land use, nitrogen fertilization, and irrigation. Intensely used 
land classes, such as villages, showed the greatest rates of increase in NDVI, more than 
twice than those of forests. These findings reveal that factors beyond climate influence 
global long-term trends in NDVI and suggest that global climate change models and 
analyses of primary productivity should incorporate land use effects. 

Keywords: NDVI; land-use; anthropogenic biomes; anthromes; global change; GIMMS3g 


1. Introduction 

Climate change research often focuses on long-term productivity trends in vegetation at 
regional [1-3] to continental and global scales [4,5]. A key variable in these studies of vegetation 
productivity is the satellite-derived Normalized Difference Vegetation Index (NDVI), that is either 
used directly as a measure of photosynthetic capacity [6], or serves as an input variable for climate 
change models predicting terrestrial net primary production (NPP [5,7,8]). On a global scale, these 
studies often ignore that NDVI may not be driven by climate alone. Because the majority of the 
Earth’s ice-free land surface has been altered by human action [9-12] and because these same areas 
account for ~90% of terrestrial net primary production [13], trends in NDVI may also be li nk ed to 
human land-use practices, such as land conversion, irrigation, and nitrogen deposition. Although 
several regional and local analyses demonstrate significant impacts on NDVI trends from changes in 
human land-use practices [14-16], global-scale analyses relating trends of NDVI to these direct 
anthropogenic effects have received little attention [5]. 

Here, we provide one of the first global analyses of direct anthropogenic effects on long-tenn trends 
in NDVI from 1981 to 2010 (however see [17] for an analysis that also include land use but on a 
different scale). First, we quantified variation in NDVI trends across “anthromes” (anthropogenic 
biomes, after [13]). Anthromes were introduced by Ellis and Ramankutty [13] and classify the earths 
land surface based on a combination of data on population density, land-use, and land cover [13]. They 
characterize land surface based on global patterns of sustained, direct human interaction with 
ecosystems. Comparisons of trends in NDVI among anthromes thus represent a useful probe of direct 
anthropogenic effects on long-tenn trends in NDVI. We also examined the effect of human-population 
density, as a surrogate for the degree of anthropogenic effects. Lastly we tested for specific 
mechanistic linkages that may underlie spatial variation in trends of NDVI among ecoregions and 
investigated effects of human land conversion (measured as percentage village, urban areas and 
croplands per ecoregion), nitrogen deposition, and irrigation. 

2. Results and Discussion 

Globally, anthromes with intensive land use (dense settlements, villages, and croplands) all had 
significant (p < 0.05) and steep increases in annual mean NDVI over a 29-year period compared to 
other anthromes. Forests exhibited marginally significant increases (J3 = 0.0006,/? = 0.06), whereas we 
found no significant trends for wildlands or rangelands (Figure 1A,C,E,G,I,K). We emphasize that 
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wildlands in the anthropogenic biomes classification used here [13] simply lump all land surfaces that 
are free or almost free of any human footprint. They thus comprise extremely varied ecosystems across 
the globe where some may show positive and other may show negative trends. Our statement only 
refers to the average across all these diverse ecosystems. 

We also calculated robust linear regression Theil-Sen (TS) slope estimators for each spatial location 
in the global NDVI data (using the GIMMS3g NDVI data). Across all anthromes, the medians of the 
per-pixel TS estimators agreed closely with the slope estimates of the annual mean (compare slope 
estimates in Figure 1 (left) with medians in Figure 1 (right)). Based on the average NDVI value in the 
respective anthromes, the TS estimators showed that NDVI of villages changed the most over the 29 
year study period (5.9% increase), followed by dense settlements (4.3%) and croplands (3.7%). With a 
2.5% increase, global net NDVI of forests increased less than half as much as villages. Global net 
NDVI increases in wildlands and rangelands were close to zero. 


Figure 1. Trends of NDVI for different groups of anthropogenic biomes (after [13]): 
Wildlands (A,B), Rangelands (C,D), Forested (E,F), Croplands (G,H), Villages (I,J) and 
Dense settlements (K,L). (Left): Annual mean NDVI and trends based on generalized least 
square models with an AR1 autocorrelation structure. P: coefficient of year, d: coefficient 
of determination, significance codes: * p < 0.05, ** p < 0.001, ns: not significant. (Right): 
Distributions of per-pixel Theil-Sen estimators, red line: median of distribution; black line 
indicates zero; als: area land surface of the globe. 
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We show that in years where global NDVI decreased, the global average NDVI of land types such 
as villages or croplands appear to have dampened minima compared to wildlands or forests (compare 
minima in Figure 1A,C,E,G,I,K). Variation induced by seasonal changes in NDVI, such as those 
associated with the interplay between human activities and monsoon dynamics [18], complicate 
interpretation of this pattern. 

In our second analysis, we tested for the effect of human-population density, as a surrogate for the 
degree of anthropogenic effects, and found strong linkages between trends in NDVI and human 
population density. Here, we used ecoregions as spatial units [19]. Ecoregions are distinct with respect 
to natural communities and resident species and thus ideal to test for spatial variation in anthropogenic 
effects on NDVI. In ecoregions with low human population densities, NDVI trends were small and 
close to zero, whereas in ecoregions with high population densities, NDVI increased significantly over 
time (Figure 2). For example, densely populated ecoregions with an average population density of 
-500 people/km 2 showed an average increase of -0.025 NDVI units over the 29 year study period (a 
6% increase compared to the global mean). Except for Australia, average annual increases in NDVI 
were also greater in more densely populated ecoregions when data were analyzed at the continental 
scale. Europe, which has the greatest average population density across ecoregions, showed the 
strongest relationship between trends in NDVI and population density (Figure 3). 
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Figure 2. Annual change of NDVI of ecoregions, based on ecoregional means of per-pixel 
Theil-Sen estimators) vs. human population density. Circle diameters are proportional to 
the size (area) of each ecoregion. Red line indicates trend based on generalized least 
squares model with exponential spatial autocorrelation structure and accounting for 
differences in sizes of ecoregions in the variance function (/?: 0.00032, p < 0.0001, A AIC: 
105, coef. determination: 0.13). For corresponding analyses by continent, see Figure 3. 
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Population density (krrf 2 ) 

Table 1 . Results from multiple regression analysis of the importance of possible drivers of 
trends in NDVI. Analyses were perfonned using ecoregional means of per-pixel TS estimators 
of trends in NDVI (n = 770, coefficient determination: 0.21, A AIC: 171, overall significance: 
< 0 . 0001 ). 



P 

/7-Value 

Intercept 

0.00017 

<0.0001 

Converted lands (% area) 

0.00040 

<0.0001 

Nitrogen (mg N/m 2 /year) 

0.00000025 

0.0001 

Irrigation (% area) 

0.00025 

0.013 


Lastly, we tested for specific mechanistic linkages that may underlie spatial variation in trends of 
NDVI among ecoregions. In particular, we tested for human land conversion (measured as percentage 
village, urban areas, and croplands per ecoregion), nitrogen deposition, and irrigation and found that all 
three of these factors had significant effects. The percentage converted lands per ecoregion, the percent 
irrigated area per ecoregion, and the amount of nitrogen deposition all were positively related to 
increasing trends of NDVI (Table 1). Together these three factors explained 21% of the variation in 
trends in NDVI. 
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Figure 3, Annual change of NDVI of ecoregions (Africa (A), Asia (B), Australia (C), 
Europe (D), North America (E), South America (F), based on ecoregional means of per-pixel 
Theil-Sen estimators) vs. log 10 of human population density. Circle diameters are 
proportional to the size (area) of each ecoregion. Red line indicates trend based on 
generalized least squares model with exponential spatial autocorrelation structure and 
accounting for differences in sizes of ecoregions in the variance function. P: coefficient of 
log 10 of human population density, d: coefficient of determination, significance codes: 

*** p < 0.0001, ns: not significant. 
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Table 2. Additional models relating trends in NDVI to (A) population density and (B) to 
nitrogen deposition, irrigation, and converted lands in order to ensure robustness of results 
against potential methodological concerns (see main text). 

A. Models Relating Trends in NDVI to Population Density 
Excluding the Moist Tropical Broadleaf Biome 

P p-Value 

Intercept -0.000054 0.25 

LoglO Population (m 2 ) 0.00036 <0.0001 

Using Standardized TS Estimators 

Coefficient of determination: 0.10, A AIC: 80, overall significance: <0.0001 

P /7-Value 

Intercept 0.000043 0.68 

LoglO Population (m 2 ) 0.00069 <0.0001 

Using Fishnet as Spatial Unit 

Coefficient of determination: 0.059, A AIC: 49, overall significance: <0.0001. 

P /7-Value 

Intercept 0.00024 <0.0001 

LoglO Population (m 2 ) 0.00023 <0.0001 

B. Models Relating Trends in NDVI to Nitrogen Deposition, Irrigation, and Converted Lands 

Excluding the Moist Tropical Broadleaf Biome 

Coefficient of determination: 0.26, A AIC: 164, overall significance: <0.0001. 


Intercept 

Converted lands (% area) 
Nitrogen (mg N/m 2 /year) 
Irrigation (% area) 

P 

0.000068 

0.00046 

0.00000023 

0.00035 

/7-Value 

0.044 

<0.0001 

0.0022 

0.0028 

Using Standardized TS Estimators 

Coefficient of determination: 0.15, A AIC: 122, overall significance: <0.0001. 


P 

/7-Value 

Intercept 

0.00029 

0.0002 

Converted lands (% area) 

0.00053 

0.024 

Nitrogen (mg N/m 2 /year) 

0.00000043 

0.0075 

Irrigation (% area) 

0.0010 

0.0001 

Using Fishnet as Spatial Unit 

Coefficient of determination: 0.11, A AIC: 88, overall significance: 

<0.0001 


P 

/7-Value 

Intercept 

0.00022 

<0.0001 

Converted lands (% area) 

0.00037 

0.0003 

Nitrogen (mg N/m 2 /year) 

0.00000015 

0.020 

Irrigation (% area) 

0.00031 

0.0064 


These results were robust to several methodological concerns and alternative hypotheses. For 
example, analyses excluding the moist tropical biome yielded similar results, suggesting the uncertainty of 
remote sensing estimates in the moist tropics [20] did not drive our findings (Table 2). Results were 
also robust to the hypothesis that humans have settled preferentially in areas with high NDVI 
(alternatively, that areas with high NDVI had greater absolute increases in NDVI) because 
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standardizing the per-pixel TS estimators by the mean NDVI did not significantly change model 
outcomes (Table 2). Finally, results were robust to the spatial unit of analyses, because substituting the 
biologically defined ecoregions with square grid cells of length 500 km also did not substantially 
change model outcomes (Table 2). 


Figure 4. (A) Trends in NDVI across the globe, 1981-2010. Ecoregional [19] extremes for 
NDVI increase (defined as the 5% of land surface with the fastest increases in NDVI, 
n = 73 ecoregions) are in red, whereas ecoregional extremes for NDVI decrease (defined as 
the 5% of land surface with the fastest decreases in NDVI, n = 38 ecoregions) are in blue; 
(B) Boxplots contrast the ecoregional extremes in A for increases (n = 73) and decreases (n 
= 38) in terms of NDVI trends, population density, percent converted lands [13], percent 
irrigated lands, and nitrogen deposition. 



Ecoregions with the fastest increases in NDVI were predominately located in temperate Europe, the 
Mediterranean, the Sudano-Sahelian belt in Africa, western India, and northern China (Figure 4A). 
Ecoregions with the fastest decreases were predominately located in boreal forests of North America, 
as well as dry forests in South America (mostly Paraguay) and drylands of northeastern Australia. 
When comparing ecoregions with the fastest NDVI increases (the 5% of land surface with the fastest 
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increases in NDVI, n = 73 ecoregions) against those showing the largest NDVI decreases (the 5% of 
land surface with the fastest decreases in NDVI, n = 38 ecoregions), we found that the ecoregions 
exhibiting the fastest NDVI increases had a greater average population density (median of 94 vs. 
7 persons per km 2 , Figure 4), greater percentage of converted lands (88% vs. 9%, Figure 4B), greater 
percentage of irrigated lands (49% vs. 4%, Figure 4B), and greater deposition of nitrogen (573 vs. 209 
mg N/m 2 /year, Figure 4B). 

Ecoregions with the fastest changes in NDVI corresponded well with decreases and increases in 
NDVI shown in other global and regional studies [5,21]. Decreasing trends in NDVI have been 
reported for ecoregions with the fastest decreases, like the Dry Chaco region in South America [14], 
northeastern Australia [22], and boreal Canada and Alaska [23,24], Likewise, increasing trends in 
NDVI have been reported for ecoregions with the fastest increases across Europe [25,26], northern 
China [27], central India [2], and the Sudano-Sahelian belt in Africa [28,29]. 

Although trends in NDVI in some of these ecoregions, especially those with negative trends in 
remote areas like boreal Canada and Alaska or northeastern Australia, have been predominately 
attributed to climatic effects [22-24], there is growing evidence for direct anthropogenic effects in 
almost all regions with significant human populations. For example, increasing trends in NDVI in 
India and northern China have been related to increasing rates of fertilization and irrigation [2,28]. 
Likewise, while previous studies have chiefly related the Sahelian greening trend to climate, e.g., [7], 
new investigations point to land conversion to croplands as one of the main causes for changes in 
NDVI [29]. In Europe, nitrogen deposition may be one of the main drivers of carbon sequestration in 
forests and — together with reforestation and climate changes — may well explain the large increases in 
NDVI observed across the continent [25,26,30]. We explain 35% of the variability in trends in NDVI 
by population density alone (Figure 3D), supporting the hypothesis that anthropogenic drivers are 
contributing to NDVI increases in Europe. There are also examples for direct anthropogenic effects in 
ecoregions with rapidly decreasing NDVI. Results of a recent study on deforestation undertaken by the 
UN-FAO [31] show that forest loss in South America accounts for more than 50% of the global forest 
area loss in the past two decades. One region particularly affected is the Dry Chaco in Paraguay, a 
region with rapid land conversion from forest to soybean plantation [32], identified by our analyses as 
an extremum for decreasing NDVI. 

3. Experimental Section 

3.1. NDVI Data and Pre-Processing 

We used NDVI data from the Advanced Very High Resolution Radiometer (AVHRR) instrument 
developed into the GIMMS3g dataset (Global Inventory Modeling and Mapping Studies) [33,34], 
Previous versions of the GIMMS data set have been used extensively for testing temporal trend in 
NDVI [35], have been tested against ground controls points [36] and are thought to be the best long 
tenn data set for analyzing long tenn trends in vegetation dynamics [37]. The current GIMMS3g 
dataset represents an improvement over previous AVHRR NDVI data sets because it was coprocessed 
with SeaWiFS and MODIS Terra data [20]. The SeaWiFS and MODIS instruments have spectral 
bands specifically designed for vegetation monitoring, better atmospheric correction algorithms, 
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reduced geometric distortions, and improved radiometric sensitivity [38]. MODIS NDVI has been 
collected over the past 13 years from the Terra platform and the past 11 years from MODIS on the 
Aqua platform. For the respective MODIS overlap periods, the per-pixel trends in NDVI from 
GIMMS3g agree well with those from Terra and Aqua MODIS, except in the moist tropics, which are 
difficult for remote sensing due to humidity and cloud cover [20], MODIS Terra and Aqua NDVI data 
do not compare well for the tropics either. Using the longer-tenn GIMMS3g dataset, we conducted the 
analyses discussed below twice, once with all terrestrial pixels and once excluding the moist tropics. 
Since our analyses used yearly means (see below) and the studies cited above did not compare yearly 
means of NDVI between MODIS and GIMMS3g, we performed additional analyses comparing annual 
means between MODIS terra and Aqua, and the GIMMS3g data set. We sampled 34 1 -degree areas 
across the globe and found a high linear correlation between the MODIS and GIMMS data sets. A 
Gaussian mixed model between annual means of MODIS Terra and GIMMS3g accounting for 
sampling area with a random effect showed very close relationships between both data sets ([> = 0.74, 
t = 24.57, N = 298, coefficient det. = 0.99) as did the same comparison with MODIS Aqua (fi = 0.75, 
t = 25.79, N = 247, coefficient det. = 0.99). If the data was standardized by the mean NDVI across 
years of each area to make areas more comparable, the relationship was still strong. MODIS Terra vs. 
GIMMS3g: (J3 = 0.68, t = 21.25, N = 298, coefficient det. = 0.65) and MODIS Aqua vs. GIMMS3g: 
( J3 = 0.69, t = 22.62, N = 247, coefficient det. = 0.65). 

The GIMMS3g dataset of bi-monthly NDVI composites provides global coverage at a spatial 
resolution of 8 km x 8 km. For our analyses, we used the 672 composites for the 29 year period 
1982-2010. We rescaled the data to the original floating-point NDVI scaling of-1 to 1, we identified 
and removed falsely included water pixels with a correction layer based on vector layers of continental 
coastline, lakes and reservoirs [39], and we applied a minimum threshold that set all pixels with NDVI 
values less than 0.05 to 0.05. NDVI values below 0.05 indicate non-vegetated areas like bare soil, water, 
snow, and ice. This 0.05 threshold is based on an analysis of stable pixels in the Sahara Desert that had 
an average NDVI of 0.094 and has been used previously in studies examining trends in NDVI [40]. 

3.2. Trend Analysis & Anthropogenic Effects 

We estimated linear trends of annual means of NDVI for the 29 point time series from 1982 to 2010 
for each spatial location in the data set (n = 1,975,140 pixels) using the Theil-Sen (TS) median 
estimator (function mblm in package mblm in R). The TS estimator represents a robust linear 
regression method that calculates the median slope among all pairs of observations. Commonly used 
for studies in vegetation trends based on NDVI time-series [20,41,42], the TS estimator is relatively 
insensitive to outliers and is accurate even for skewed and heteroscedastic data. 

To relate these estimates to anthropogenic effects we perfonned three different analyses. First, we 
evaluated trends in NDVI for each of the six groups of anthropogenic biomes identified by Ellis and 
Ramankutti [13]; wildlands, rangelands, forests, croplands, villages, and dense settlements. We calculated 
histograms and medians for the TS estimators for each group of anthropogenic biomes and the global 
trend for each group based on the yearly average of NDVI for all pixels of a group using generalized 
least squares (GLS). We incorporated the temporal correlation structure with a 1st order autoregressive 
moving average. For all GLS analyses, we calculated the coefficient of detennination based on the 
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improvement of the log-likelihood from the fitted null model (intercept only) to the full fitted model 
(function r.squaredLR in R package MuMIn, [43]). 

Our second and third analyses were carried out on the scale of ecoregions. We used the World 
Wildlife Fund’s classification scheme, which distinguishes ecoregions [19] that are distinct with 
respect to natural communities and resident species. We first tested whether human-population density, 
as a surrogate for the degree of anthropogenic effects, was related to NDVI trends. We estimated 
human population density for each ecoregion as the mean of the 1990 and 2010 population estimates [44], 
We next tested for the effects of land conversion, nitrogen deposition, and irrigation on NDVI trends 
by performing multiple regressions that included these factors as predictors. As a proxy for converted 
lands, we used the percent area covered by dense settlements, villages, and croplands per ecoregion [13]. 
For nitrogen deposition we used the mean nitrogen deposition per ecoregion in 1993 ([45], an average 
from the 1980s to the 2010s would have been preferable but was not available). For irrigation we 
calculated the percentage of irrigated area per ecoregion based on the UN-FAO irrigation map [46]. 

For both analyses we used the ecoregional means of the per-pixel TS estimators as the response 
variable. We used GLS accounting for spatial autocorrelation with an exponential spatial correlation 
function and accounting for differences in sizes of ecoregions using the inverse of ecoregion size as a 
weight in the variance function. For the multiple regression, we tested predictors for collinearity and 
used a cutoff correlation coefficient of 0.7 [47]. 

For robustness, we first repeated the analyses excluding all pixels located in the moist tropical 
broadleaf biome, because of unresolved and ongoing debates about the validity of remote sensing 
estimates in moist tropical regions [20,48] and the poor agreement between (1) GIMMS3g and MODIS 
Terra NDVI and (2) MODIS Terra and Aqua NDVI data in the moist tropics found by Fensholt and 
Proud [20]. Second, to ensure that our results do not stem from humans choosing areas with greater 
NDVI or from areas of greater NDVI exhibiting greater absolute increases of NDVI, we repeated our 
analyses by standardizing the per-pixel TS estimators by the mean NDVI value of each pixel. Third, 
even though ecoregions represent ecologically distinct units and as such are appropriate geographic 
units for examining the effects of covariates on trends in NDVI, we also repeated our analyses 
substituting squared grid cells of 500 km by 500 km for ecoregions. Lastly, to examine whether our 
results hold true also within continents and not just on the global scale, we perfonned the analyses for 
population density separately for each continent. Because some continents have only relatively few 
ecoregions, we did not perform the by-continent analyses for the multiple regression analyses. 

4. Conclusions 

In conclusion, our study provides one of the first global synopses quantifying the effect of the 
Earth’s human footprint on trends in photosynthetic capacity. Globally, more than 20% of the 
variability in NDVI trends can be associated with human land-use practices. Consequently, global 
climate change models and analyses of primary productivity should incorporate direct anthropogenic 
effects. Land use practices in densely populated areas, and other areas featuring major human impacts 
through land conversion, irrigation, and nitrogen deposition, may contribute more to global changes in 
photosynthetic capacity than previously anticipated. 
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